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We have calculated complete spectra of the staggered Dirac operator on the lattice in quenched 
SU(3) gauge theory for l3 — 5.4 and various lattice sizes. The microscopic spectral density, the 
distribution of the smallest eigenvalue, and the two-point spectral correlation function are analyzed. 
We find the expected agreement of the lattice data with universal predictions of the chiral unitary 
ensemble of random matrix theory up to a certain energy scale, the Thouless energy. The deviations 
from the universal predictions are determined using the disconnected scalar susceptibility. We find 
that the Thouless energy scales with the lattice size as expected from theoretical arguments making 
use of the Gell-Mann-Oakes-Renner relation. 
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The low-lying eigenvalues of the Dirac operator are of 
great importance for the understanding of spontaneous 
chiral symmetry breaking in an infinite volume Q. On 
the lattice, however, one is always working at finite vol- 
ume. Therefore, it is important to know how the thermo- 
dynamic limit is approached. It was shown by Leutwyler 
and Smilga S that in the domain 



1/A < L < 1/m^ 



(1) 



where A is a typical hadronic scale, L is the linear extent 
of the Euclidean box, and is the pion mass, the low- 
energy behavior of QCD can be described by a simple 
effective partition function whose existence imposes cer- 
tain constraints on the eigenvalues of the Dirac operator. 
The spectrum of the Dirac operator in the domain (|l|) 
has been successfully predicted by chiral random matrix 
theory (RMT) The only ingredients of the calcu- 

lation are the global symmetries of the theory and the 
assumption that chiral symmetry is spontaneously bro- 
ken. It has recently been shown that there is an overlap 
between the domain of validity of chiral RMT and of chi- 
ral perturbation theory, and that in this overlap region 
the two approaches yield the same results Q . 

The results so obtained provide analytical information 
on the way in which the thermodynamic limit is ap- 
proached. They are universal in the sense that they do 
not depend on the precise values of the parameters of the 
theory, i.e., of the simulation parameters on the lattice. 
However, the domain of validity of the universal results 
does depend on the parameters. The energy scale Armt 
up to which RMT applies, i.e., the Thouless energy, fol- 
lows from the upper bound on L in relation (|^) and the 
Gell-Mann-Oakes-Renner relation, m'^f^ = 2mE, where 

= F^/\/2 is the pion decay constant, E is the absolute 
value of the chiral condensate (tpip), and m is a valence 
quark mass. It is thus determined by 



2t2 



Armt/A cx f^L 



(2) 



where A = l/p{0) = tt/{VJ^) is the level spacing at 
zero. Here, V ^ L'^ denotes the four-volume, and 
pW — {J2n^i^ ~ ^n))A is the spectral density of the 
Dirac operator averaged over gauge field configurations 
A. The relation between p{Q) and S is given by the 
Banks-Casher formula, 7rp(0) = VII [|l|. 

The aim of this paper is (i) to test the universal pre- 
dictions of chiral RMT for the distribution and correla- 
tions of the low-lying Dirac eigenvalues and (ii) to check 
the prediction of Eq. (||) for the Thouless energy, using 
lattice data computed in quenched SU(3) gauge theory 
with the staggered Dirac operator. Point (i) has previ- 
ously been considered in quenched SU(2) |P,|lOt, in SU(2) 
with dynamical fermions pT| ], in quenched SU(3) in three 
dimensions [ p^ , in U(l) in two dimensions p^ , and, 
very recently, in quenched SU(3) in four dimensions [Q. 
Point (ii) has previously been tested in quenched SU(2) 
|]l5| . All these investigations were done with the stag- 
gered Dirac operator except for Ref. jl^] in which the 
fixed point Dirac operator with respect to a renormaliza- 
tion group transformation was used. Since real QCD has 
three colors, SU(3) in four dimensions is clearly the most 
important case. 

The Euclidean Dirac operator in the continuum is 
given by ip = ip+gt"'j^'^ , where the t° are the generators 
of the gauge group. The operator ip is hermitean with 
real eigenvalues. It anticommutes with 75 and, therefore, 
all nonzero eigenvalues come in pairs ±A„ with eigen- 
vectors V'liiTsV'ii- There can also be zero modes which 
are either left-handed or right-handed. The topological 
charge of a given gauge field configuration is equal to the 
difference in the number of left-handed and right-handed 
zero modes. On the lattice, the staggered Dirac operator 
reads 



1 



(3) 



where U and 77 denote the hnk variables and the staggered 
phases, respectively. 

The claim is that the distribution and the correlations 
of the small eigenvalues of ip are described by universal 
fimctions which can be computed, e.g., in chiral RMT 
. The distribution of the low- lying eigenvalues is en- 
coded in the spectral one-point function near zero virtu- 
ality, the so-called microscopic spectral density defined 

by i 



Ps{z) = lim jj^ p I — 1 



(4) 



Similarly, one considers the microscopic limit of the two- 
point cluster function. 



with 



r2(Ai, A2) = p(Ai)p(A2) - i?2(Ai, A2) , 



(5) 



(6) 



where i?2(Ai,A2) is the two-point spectral correlation 
function, i.e., the probability density that one eigenvalue 
is at Ai and another at A2, all other eigenvalues being 
unobserved. 

For SU(3) with the staggered Dirac operator, the rel- 
evant symmetry class in the framework of chiral RMT 
is the chiral unitary ensemble |16|. In the following, 
we briefly summarize analytical results for this ensem- 
ble which are of relevance for the present work. The 
microscopic spectral density is given by H 



(7) 



with the Bessel function J and fi — Nj + where 
Nf and i> denote the number of massless flavors and the 
topological charge, respectively. The distribution of the 
smallest eigenvalue for Nf — v = reads [ pTf 



-P(Amin) — 



A. 



(8) 



The two-point cluster function in the microscopic limit 
is given by 



T2{Z1, Z2) 



Z1Z2 



ZlJf_i+lizi)J^,{z2) - Z2Jfj.izi)J^+i{z2) 



Z2 



(9) 



The quantities in Eqs. <Q) through (|9[) do not contain any 
free parameters. For a comparison with lattice data, the 
energy scale is determined by the parameter VT, which 



TABLE I. Simulation parameters and extrapolated value 
of VT, for 13 = 5.4. 



L 


conf. 


VT, 


4 


35337 


225 ± 7 


6 


11748 


1207 ± 28 


8 


2635 


3918 ± 58 


10 


1059 


9429 ± 155 



is obtained from the data by extracting p(0) and apply- 
ing the Banks-Casher relation, 7rp(0) = VJ^. Thus, the 
comparison between lattice data and the predictions of 
Eqs. (j^) through (||) is parameter-free. (Strictly speak- 
ing, on finite lattices a spontanous breaking of chiral sym- 
metry cannot occur and p(0) is zero. The latter quantity 
must, therefore, be determined by extrapolating p(A) to 
A = 0. In practice, this extrapolation presents no difh- 
culties.) 

We now turn to the details of our numerical simula- 
tions. They were done in quenched SU(3) gauge the- 
ory with (3 — G/g^ = 5.4 on lattices of size V = L'^ 
with L = 4, 6, 8, 10. The boundary conditions are peri- 
odic for the gauge fields and periodic in space and anti- 
periodic in Euclidean time for the Dirac operator. The 
gauge field configurations were generated using a com- 
bined Metropolis and overrelaxation algorithm on the 
link variables. Two consecutive configurations are sep- 
arated by at least 30 runs of one Metropolis sweep with 
three hits and 20 overrelaxation sweeps using Creutz's 
method jl^. The complete spectrum of the staggered 
Dirac operator was then calculated using the CuUum- 
Willoughby version of the Lanczos algorithm for the ma- 
trix of —p^. This operator couples only even to even and 
odd to odd lattice sites. Both blocks have the same eigen- 
values. Hence it is sufficient to consider only even lattice 
sites. The eigenvalues of —p"^ were checked against the 
identity Tr(— = 3^ which was fulfilled with relative 
accuracy 10~®. The total number of diagonalizcd config- 
urations and the extrapolated values of Trp{0) = VS are 
shown in Table |. 

The lattice data for Ps{z) and P(Ainin) are compared 
with the predictions of Eqs. (|^) and (^) in Fig. |l|. In 
Eq. (|7|), we have used ^ = 0. Clearly, Nf — since we 
consider the quenched approximation. The fact that — 
is less obvious. The prediction of Eq. (Q) is restricted 
to sectors with definite topological charge. Therefore, 
one should compute the topological charge of each gauge 
field configuration and compare the lattice data in each 
topological sector with the prediction of Eq. (|^). How- 
ever, Eq. ^ assumes that for 1/ ^ the Dirac operator 
has exact zero modes. This is not the case for staggered 
fermions on the lattice where at finite lattice spacing a 
the would-be zero modes are shifted by an amount pro- 
portional to jl^]. For the value of (3 we used, a is still 
relatively large so that no zero modes are present. This 
explains why the lattice data are consistent with Eq. (Q) 



2 




2 4 6 8 10 12 14 




2 4 6 



10 12 14 



2 4 6 8 10 12 14 




2 4 6 8 10 12 14 






0.0 0.0005 ^ 0.001 



2,500 




0.0002 0.0004 



FIG. 1. Microscopic spectral density (left) and distribution 
of the smallest eigenvalue (right) of the staggered Dirac oper- 
ator in quenched SU(3) for /3 = 5.4 and V = L"^ with L = 4, 6, 
8, 10. The histograms represent the lattice data, the dashed 
curves are the predictions of Eqs. (^ (with = 0) and (^). 
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FIG. 2. Microscopic two-point cluster function T2{zi,Z2) 
with Z2 at the second maximum of ps{z) for L — 4 and L — 6, 
respectively. The histograms represent the lattice data, the 
dashed curves are the prediction of Eq. (^) for /i = 0. 

for = 0, as seen in Fig. ^ Very similar results for 
different /3 were very recently presented in [T^ ] . 

The agreement between the lattice data and the uni- 
versal predictions is quite satisfactory, also for the two- 
point cluster function in the microscopic limit which we 
have plotted in Fig. |^ along with the prediction of Eq. (|^) 
for fi ^ Q. A fixed value of Z2 = 4.68 (corresponding to 
the location of the second maximum of Ps{z)) was chosen. 

The quantity T2 (2:1,22) is interesting since it enters 
in the calculation of the disconnected scalar susceptility 
which, in turn, can be used to determine the Thouless 
energy, i.e., the scale Armt above which the lattice data 
deviate from the universal predictions of Eqs. ^ through 
(^. In terms of the Dirac eigenvalues, this quantity is 
defined as I20II 



disc 
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^ {i\k + m){iXe + m) 
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1 



N 

y 



(10) 



where V is the number of lattice sites, N the number 
of eigenvalues and m a valence quark mass, respectively. 
The average is over independent gauge field configura- 
tions. Eq. (|l^) can be rewritten in terms of integrals in- 
volving the spectral one- and two-point functions of the 
Dirac operator. Rescaling by \/{VT?) and changing 
from TO to M = mVY,, we have 



dx 



dxdy ■ 



10 JO (x2+u2)(j/2 
x[^'(^)--^A'+l(w)Vl(")] ' 



(11) 



where in going from the first to the second line we have 
inserted the RMT results for ps and T2. The functions 
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FIG. 3. The ratio defined in Eq. (|13|) for (3 = 5.4 and vari- 
ous lattice sizes. 



FIG. 4. The ratio of Eq. (hj) plotted versus u/L'^ 



I and K are modified Bessel functions. In the case of 
^ = 0, Eq. (|l^) simplifies to 

X^'^Hu) = u'[K^Au) ~ Xo'(")] [loi^) - Ifi^)] ■ (12) 

In order to compare the lattice results for 1^'^""= ob- 
tained from Eq. with the RMT prediction of Eq. ^) 
we introduce the variable llSfl 



ratio = 



^disc __disc 

Alatticc ARMT 

^disc 

ARMT 



(13) 



which is plotted in Fig. ^. This ratio should be close to 
zero in the domain of validity of the RMT predictions 
and deviate from zero at some value of u = mrmt which 
corresponds to the Thouless energy. (The deviations of 
the ratio from zero for very small values of u are artefacts 
of the finite lattice size and of finite statistics. This point 
was discussed in Ref. [jTsly 

The prediction of EqT"^) is that Armt / A should scale 
with L^. If we express Armt in terms of urmt = 
ArmtV^S — ttArmt/A, urmt should also scale with L^. 
To check this predicted scaling behavior we have plotted 
the ratio of Eq. ( p^ ) as a function of u/L'^ in Fig. ^. We 
observe that all data fall on the same curve, confirming 
the prediction of Eq. (^ with regard to the scaling with 
L^. Since we have only considered one value of (3, we 



cannot check the scaling with From Fig. |j we can 
read off Armt/ A ~ 0.04i^ in lattice units. 

In summary, we have shown that the distribution and 
the correlations of the low-lying eigenvalues of the stag- 
gered Dirac operator in quenched SU(3) are described 
by universal functions up to a certain energy scale, the 
Thouless energy. The latter quantity was determined 
using the disconnected scalar susceptibility, and the pre- 
dicted scaling with was confirmed. It would be of 
great interest to extend the present study to dynamical 
fermions for which analytical results are also available 
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